Controlling crystal self-assembly using a real-time feedback scheme 
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We simulate crystallisation of hard spheres with short-ranged attractive potentials, as a model 
self-assembling system. We show how measurements of correlation and response functions during 
assembly can be used to tune the interaction parameters as assembly proceeds, in order to obtain 
high-quality crystals. The method we use is independent of details of the interaction potential 
and of the structure of the final crystal - we propose that it could be applied to a wide range of 
self-assembling systems. 



In self-assembly [TJ [5] , simple components like colloids 
and biomolecules come together spontaneously, forming 
ordered structures such as viral capsids [3H], crystals [5]- 
[9], or DNA origami [10]. Recent developments in self- 
assembly include the experimental synthesis of particles 
with specific controllable interactions [71I51 ITTHT3"] . as well 
as theoretical and computational demonstrations of the 
range of ordered phases that such particles can form 
fT7] , However, even when ordered states are stable, aggre- 
gation into disordered states often frustrates the dynam- 
ical self-assembly process. As a result, effective assem- 
bly typically involves microscopic reversibility: if bonds 
are both made and broken during self-assembly then de- 
fects are annealed naturally, producing an ordered final 
state O HI [9l IT8H20"] . The twin requirements of a stable 
ordered structure and the reversibility of bonding usu- 
ally mean that assembly is effective only when interac- 
tion parameters are tuned to lie within a narrow 'optimal' 
range [3 E USHU] • 

If interactions between particles can be manipulated 
during the assembly process, this may facilitate the op- 
timization of assembly conditions. For example, a slow 
cooling process is used in hierarchical assembly of DNA 
origami |10) . and light exposure in repeated pulses may 
allow control of nanoparticle aggregation [2H [53] . How- 
ever, introducing time-dependent interactions in self- 
assembly comes at a high price in complexity, since there 
are very many possible protocols to be considered. Here, 
we develop a scheme where we vary the strength of 
interactions between particles with time, according to 
measurements of the reversibility of particle bonding, 
as assembly proceeds. This cycle of measurements and 
changes in interaction strength forms an elementary feed- 
back loop. Our results provide a proof of concept for this 
approach, which represents a systematic way of exploring 
the wide range of possible assembly protocols. 

We consider crystallisation of hard spheres with short- 
ranged attractive interactions as a model self-assembly 
process. Such models have been used to represent col- 
loidal particles interacting through depletion forces [9] 
H9] [24] , and as simplified descriptions of protein crystalli- 
sation 25J. Using computer simulations, we show that 



(i) our feedback scheme quickly locates the regime where 
assembly is effective, and (ii) it forms higher-quality crys- 
tals than assembly with fixed (time-independent) inter- 
actions. We emphasise that the feedback scheme does 
not rely in any way on the structure of the ordered (crys- 
tal) state - we expect that it can be applied to other 
self-assembly processes with minimal modification. 

Model - We consider N = 10, 000 spherical parti- 
cles with hard cores of diameter a, at volume fraction 
4> = 0.04, interacting through a pair potential U(r) which 
is a square well of depth e^, and range £ = O.Ict. The tem- 
perature is T and the particles evolve by overdamped 
Langevin dynamics with (bare) diffusion constant D; we 
define a time unit tq = a 2 /D, of the order of a Brow- 
nian time. The dynamical evolution is implemented by 
single particle Monte Carlo (MC) moves with parameters 
as in [5]: hence To ~ 270 MC sweeps. 

For simulations with fixed (time-independent) interac- 
tion parameters and times up to 10 To, we find appre- 
ciable crystallisation only for 2.3 < e^jT < 3, with the 
highest yield obtained towards the centre of this range 
(see also below). We use 'No- Feed' to indicate these re- 
sults, while 'With-Feed' indicates results obtained using 
the feedback scheme. 

Measurements - The feedback scheme is based on esti- 
mates of the reversibility of particle bonding, using corre- 
lation and response functions, out of equilibrium [26ll30j . 
We have shown previously that these measurements are 
sensitive to the reversibility of assembly processes, and 
are therefore strongly correlated with the effectiveness of 
self-assembly (9][27J[3T]. To make use of this correlation, 
we perturb the dynamics of the system, writing the en- 
ergy as E = \ T, p 0-+h v ) E p where E p = E P ' U{\r p ~r p ,\) 
is the energy of particle p and h p is a perturbing field ap- 
plied to that particle. We define correlation and response 
functions, 

C(Mw) = (E p (t)E p (t w )) - (E p (t))(E p {t w )) (1) 
x{t,t w ) = ^^{E p {t)) (2) 

In practice, it is convenient to estimate x by taking 
h p = +h for half of the particles in the system and 
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FIG. 1: Operation of the feedback loop on short time scales 
for three 'With-Feed' simulations with different initial bond 
strengths eJJ" . (Top) Time-dependence of the bond strength 
£b- Measurement and relaxation stages are indicated for the 
run with e™ 1 ' JT = 7.0. The shaded region indicates where 
'No-Feed' simulations found appreciable crystallisation. The 
inset shows the time-dependence of the perturbation strength 
h for the run with /T = 7.0. (Bottom) Average number 
of bonds per particle (n p (t)) = fiff . 



hp = —h for the other half. To accurately estimate the 



linear response defined in ([2|), we take h = hT/eh with 
h = 0.1. At equilibrium, correlation and response func- 
tions are related by the fluctuation-dissipation theorem 
(FDT), x cqm (Mw) = C cc i m (i,i)-C cc i m (M w ). Of course, 
self-assembling systems are far from equilibrium so FDT 
does not hold: the deviations from FDT behaviour indi- 
cate the balance between the microscopically reversible 
particle bonding and the macroscopically irreversible self- 
assembly [31] . 

Overview of feedback scheme - The feedback scheme 
consists of alternating 'measurement' and 'relaxation' 
stages. In measurement stages, we perturb the parti- 
cles with fields h p and we measure correlation and re- 
sponse functions. The duration of the ith measurement 
stage is Aj, which depends on the behaviour of correla- 
tion and response functions during that stage. At the 
end of each measurement stage, the bond strength e b is 
updated, according to the measurements that have been 
made. Each measurement stage is followed by a relax- 
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FIG. 2: Behaviour of the correlation and response func- 
tions for a 'With-Feed' simulation with e™ /T = 7.0 (see also 
Fig.Q. 



ation stage of the same duration, Aj, during which h = 0. 
At the end of each relaxation stage, the next measure- 
ment stage begins. The scheme begins with a relaxation 
stage of duration A = IOto- The behaviour during the 
first few stages of the feedback loop is summarised in 
Fig. [T] The first key result of this paper is that for a 
range of initial bond strengths, the system quickly con- 
verges to the (narrow) range of £b/T where assembly is 
effective. (The time scales required for effective crystalli- 
sation are much longer than those shown here: typically 
at least 5, 000 - 10, 000t o .) 

Duration of measurement stages - Measurements of 
reversibility of bonding during assembly must be made 
on time scales where significant bond-making and bond- 
breaking has taken place. These time scales are not 
known a priori: they depend on the particle interactions 
and they also change significantly during the assembly 
process. We therefore determine these time scales adap- 
tively, using information contained in the correlation- 
response measurements themselves. If the ith measure- 
ment stage begins at time ti, it is convenient to nor- 
malise correlation and response functions as C(t,ti) = 
C(t,ti)/C(t,t) and x(t,U) = x{t,U)/C{t,t). We also 
define X(t,t % ) = x(t,k)/[C(t,t) - C(t,ti)]: for a re- 
versible (equilibrated) system then X — 1; for irreversible 
aggregation we typically find X rj 0. The ith mea- 
surement stage ends when t is large enough that either 
C(t,ti) < C* or X(t,ti) < X*, where C* = 0.3 and 
X* = 0.6 are parameters associated with the feedback 
scheme. 

To illustrate this, consider Fig. [2j where the behaviour 
of C and x are shown for one of the 'With-Feed' simula- 
tions from Fig. [lj The functions 1 — C(t, ti) and x(t, *w) 
would coincide at equilibrium: their ratio is equal to X. 
(Note C(t,t) = 1 by definition.) In the first measure- 
ment stage, the bonds between particles are strong and 
the bonding almost irreversible - thus, x is small and the 
ratio X is much less than its cutoff X* . As a result, the 
measurement stage is quickly terminated and the bond 
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FIG. 3: Configurations from 'No-feed' (left) and 'With-feed' (right) simulations at t = 10 4 ro, and averaged measures of 
crystallinity at this time. In the snapshots, particles in fcc/hcp environments are coloured orange /purple. Particles in other 
environments are rendered at half their actual diameter for visual clarity. The data in the table are averaged either over 8 
independent simulations ('No-Feed') or over the M = 8 systems associated with single 'With-Feed' simulations. 



strength reduced (recall Fig. [T]) . In the second and third 
measurement stages, the responses are larger, indicating 
that bonding is more reversible. These stages terminate 
when C has fallen to C* , indicating that the system has 
decorrelated significantly from its state at time tf. by fix- 
ing a value for C* we ensure that different measurements 
of reversibility can be compared on an equal footing even 
when the absolute time scales for bonding and unbonding 
are changing. 

Update rule for e D - At the end of each measurement 
stage, we update the bond strength e b by a dimensionless 
multiplicative factor that depends on X. For simplicity, 
we take 



new _old 



l + -(X-X a ) 

7t 



(3) 



Here Xq is the target for the response strength X. We 
take Xq = 0.8, indicating that for ideal assembly devia- 
tions from reversible behaviour should be significant but 
not too large: this has been shown to correlate with good 
assembly in several self-assembling systems [91 |571 I3"T] . 
Also, a — 1 determines the sensitivity of the feedback 
loop, and tj is a damping parameter. We take 71 = 1 
for the first measurement phase; at the end of subse- 
quent measurement phases then 7.; is increased by 1 if 
the response X is within a tolerance 5X of Xq. We take 
5X = 0.1. 

This feedback scheme does require several parameters, 
such as the target Xq, cutoffs (C* ,X* ,SX) and sensitiv- 
ity a. The values we have chosen for these parameters re- 
flect our understanding of the physics of reversible bond- 
ing and unbonding in self-assembly [HI HH1 123 GH] ■ We 
emphasise that these physical considerations are not spe- 
cific to the model system considered here, nor indeed to 



crystallisation: we anticipate that operating this scheme 
with unchanged parameters on different self-assembling 
systems will give similarly effective results. Also, we do 
not expect the behaviour shown here to be sensitive to 
small changes in the feedback parameters. 

Implementation of feedback scheme - In general, accu- 
rate estimation of x requires good statistics, and our im- 
plementation of the feedback loop is designed to achieve 
this. At the beginning of each measurement stage, the 
configuration of the system is saved and the measure- 
ment stage is simulated twice, once with randomly chosen 
hp and the second time with fields h p that are opposite 
to the initial choice. The response is averaged over the 
two simulations, minimising random errors arising from 
the specific choice of the h p . In addition, we simulate 
M — 8 systems in parallel, each with N = 10, 000 parti- 
cles. Each system has the same (time-dependent) bond 
strength e D - The averages in (JlJ and ^ run over all 
particles in all systems, and over both simulations of the 
measurement stage: this information is all combined to 
determine the duration of the measurement stage and the 
update in the value of e^. We also note that if we mea- 
sure X > 1 then we take X = 1 in ([3]): this minimises 
the sensitivity of the system to random errors in our es- 
timate of x when bonds are weak. We emphasise that 
while accurate estimation of x does require some com- 
putational overhead, the usual alternative is to perform 
many No- feed simulations over a range of £t>, in order 
to locate the good-assembly regime: that approach also 
requires significant CPU time. 

Evaluation of crystal quality - We now analyse the 
crystals that are assembled within No-Feed and With- 
Feed protocols on the relatively long time scale t = 10 4 To: 
see Fig. [3] Results at longer times are qualitatively sim- 
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FIG. 4: Behaviour of bond strength ey, at long times, for the 
'With-Feed' simulations of Fig. [IJ The dashed line shows the 
'No-feed' protocol that gave the best yield. 



ilar, although there is some slow growth in crystalline 
order as coarsening (Ostwald ripening) takes place. To 
assess the effectiveness of crystallisation, we measure lo- 
cal crystalline order using a common neighbour analy- 
sis (CNA) [3^. The CNA assigns a 4-digit signature 
to each bonded pair of particles in the system. Bonds 
with '1421' or '1422' signatures are characteristic of close- 
packed crystals: we measure the mean number of such 
bonds per particle n^, as in [S]. We also use the CNA 
to identify particles whose local environment is consistent 
with face-centred cubic (fee) or hexagonal close-packed 
(hep) order. (Fee particles have 12 bonds with '1422' 
signatures while hep particles have 6 bonds with '1421' 
signatures and 6 with '1422'.) It may be seen from Fig. [3] 
that the 'With-feed' scheme has produced a close-packed 
crystal, albeit with random stacking of fec/hep planes. 
On the other hand, the 'No-Feed' simulations give some 
clusters where five-fold packing defects are apparent, pre- 
sumably due to growth around a critical nucleus that 
lacks the symmetry of the crystal. 

We also measure the typical size of crystalline do- 
mains in the system using bond-order parameters. As 
in [531 [35], let qe (p) consist of the projection of the 
bonds of particle i onto the spherical harmonics with 
I = 6, normalised so that q%(p)* ■ %(p) — 1- Then, if 
Q = J2 P =i <1&{p) then Mq = A _1 (Q* • Q) is an estimate 
of the typical crystalline domain size. (To sec this, we as- 
sume that q e (p)* ■ qe(p') ~ 1 if particles p and p' belong 
to the same domain, with qe(p)* ■ qe(p') ~ otherwise.) 

We compare results from different protocols in Fig. [3j 
the 'With-feed' protocol results in slightly fewer crys- 
talline bonds (measured by 71142), but with significantly 
larger crystalline domains (measured by Mq). We con- 
centrate on the formation of crystals with order on large 
length scales, since applications in photonics or X-ray 
diffraction both require well-developed Bragg resonances, 
so we take A/q as our main figure-of-merit for crystalli- 



sation. By this measure, the With-feed protocol clearly 
outperforms the No-feed simulations. This is our second 
key result. We also note that there are differences in the 
crystalline yield of different 'With-Feed' runs: this comes 
partly from the choice of initial bond strength but there 
is also a stochastic component due to noise in the esti- 
mate of correlation and response functions: recall Fig. [2] 
However, the trends that we discuss here are all robust 
to these differences. 

To understand the effectiveness of the With-Feed pro- 
tocol, consider Fig. [4j On short time scales, the feedback 
scheme leads to stronger bonds than the best No-Feed 
protocol, but on longer time scales the feedback mech- 
anism reduces the bond strength. Inspection of the as- 
sembly trajectories shows that for these bond strengths, 
the short-time physics (t < IOOto) is dominated by nu- 
cleation, while the behaviour on longer times is that of 
Ostwald ripening, with smaller clusters shrinking and 
larger ones growing. In this case, the With-Feed protocol 
leads naturally to a scheme where the bonds are initially 
strong, promoting nucleation; on later time scales the 
bonds are weaker, which promotes rapid Ostwald ripen- 
ing, (this is a thermally activated process in these sys- 
tems, proceeding faster when bonds are weaker). 

This long-time behaviour within the feedback process 
arises due to properties of correlation and response func- 
tions during coarsening [26 . Briefly, the behaviour of 
crystalline clusters in this regime depends on their macro- 
scopic properties (their curvature and surface tension), so 
applying the perturbation h p to particle p does not affect 
whether the cluster containing that particle should grow 
or shrink. As a result, the response x is small in the 
long-time (coarsening) regime, and the With-Feed proto- 
col pushes the system to increasingly weak bonds. For 
very long times, it is possible that the With-Feed proto- 
col may reduce £b so far that the bonds become too weak 
and the clusters disintegrate. We did not observe this, 
although we do find a similar effect when the damping 
term 7, is not included in Certainly these possibili- 
ties need to be considered when applying this scheme to 
other systems. 

To conclude, we have introduced an algorithmic 
method for tuning interaction parameters during self- 
assembly. We emphasise the high quality of the crystals 
that it produces, and that it can be generalised straight- 
forwardly to computer simulations of other assembling 
systems. Furthermore, there is potential for applica- 
tions in experiments: attractions between colloidal par- 
ticles can now be manipulated in real-time [531 I54H59] , 
and correlation-response measurements have also been 
made [J0H32]- Still, combining these into an effective 
feedback scheme seems quite challenging. 
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